One class of integrals evaluation in magnet 

solitons theory 

Zhmudsky A. A. 
February 2, 2008 

Abstract 

An analytical-numeric calculation method of extremely compli- 
cated integrals is presented. These integrals appear often in magnet 
soliton theory. 

The appropriate analytical continuation and a corresponding inte- 
gration contour allow to reduce the calculation of wide class of inte- 
grals to a numeric search of integrand denominator roots (in a com- 
plex plane) and a subsequent residue calculations. The acceleration 
of series convergence of residue sum allows to reach the high relative 
accuracy limited only by roundoff error in case when 10 -f- 15 terms 
are taken into account. 

The circumscribed algorithm is realized in the C program and 
tested on the example allowing analytical solution. The program was 
also used to calculate some typical integrals that can not be expressed 
through elementary functions. In this case the control of calculation 
accuracy was made by means of one-dimensional numerical integration 
procedure. 

1 Introduction 

Nonlinear excitations (topological solitons 0) play an important role in the 
physics of low-dimensional magnets [|2], |3| . They contribute greatly to a heat 
capacity, susceptibility, scattering cross-section and other physical charac- 
teristics of magnets. In particular, for two-dimensional (2D) magnets with 
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discrete degeneracy it is important to take into account the localized stable 
(with quite- long life time) 2D solitons 0, [|. According to the experiments 
||], these solitons determine the relaxation of magnetic disturbance and can 
produce peaks in the response functions. 

The traditional model describes the magnet state in terms of the unit 
magnetization vector m, fh 2 = 1 with the energy function in the form 

W = J d 2 x{A(Vm) 2 + w (m)} , (1) 

where A is the nonuniform exchange constant and wo(m) - the anisotropy 
energy. 

In an anisotropic case the solution is mult i- dimensional, that is why the 
soliton structure is determined by the system of equations in partial deriva- 
tives. There are no general methods of finding the localized solutions to such 
equations and analyzing the stability of the solution. For this reason a direct 
variational method is often used. Consequently, a choice of a trial function 
plays a key role in such analysis [|, |E|, [7J. In most quite successful 

choice of trial function is as follows: 

R ( T \ 

tg- = —exp( — ) (1 + CiCOs2x), V 5 = X + C 2 sin 2x + <p , (2) 
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where r, \ are the polar coordinates in a magnetic plane, R, a, C\, C 2 and 
ifo - variated parameters. 

In papers |5|, |6|, 0] the iterative Newton method of solving the non-linear 
system of algebraic equations § was used to find the variated parameters 
providing an energy minimum. The numerical algorithm mentioned above 
results in the necessity of multiple calculation of two-dimensional integrals 
(by the angle x an d radius r) . Experience reveals that the calculation process 
could be essentially accelerated (and the calculation precision improved) if 
the analytical-numeric procedure for the integrals by radius are made before- 
hand. The typical form of these integrals is: 

f(r) exp(mr)dr 

(3) 
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[(r exp(r)) 2 + a 

where n,m - integer numbers (m < 2n), f(r) - polynomial of degree k in 
variable r. It is convenient to denote g(r) = f(r) exp(mr). Possible expres- 
sions for g(r) are rexp(r), (1 + r)exp(r), r 2 exp(r) and so on. Note, that a 
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is essential parameter - no substitution exists to eliminate it. Below we shall 
consider the case n = 1. From the one hand it is, practically, always possible 
to reduce the calculation at n > 1 to a sum of integrals at n = 1. From the 
other hand it is not difficult to modify evaluation scheme if n > 1. 

Further we show that the indicated type of integrals may be calculated 
both analytically and numerically (by the use of joint analytical and numer- 
ical methods). To obtain the desired accuracy one should find integrand 
residues (numerically) and build approximate expression (see below). 

Main analytical expressions that allow to find integrals like @) are pre- 
sented in Analytical Formulae and the C version of corresponding program 
is given in Program Realization. 

The theory of function of a complex variable gives the powerful method of 
definite integrals calculation. In the considered case the algorithm realization 
is quite problematical because the expressions for the roots in a complex plane 
could not be written analytically. Correspondingly, one can not write down 
and analytically summarize expressions for residues. 

Present paper pays attention to the fact that numerical methods (roots 
search in complex plane and calculation of residue sum) together with analyt- 
ical ones allow to use the theory of complex variable with the same efficiency 
like in the traditional consideration. 



2 Analytical formulae 

The usual way of integral @ evaluation is to continue analytically the in- 
tegrand function in some complex domain D. In this case the evaluating 
integral will be the part of the integral over the closed contour C in a com- 
plex plane. This contour comprises the real axis interval and the arc S closing 
the integration contour. So, the solution of the problem can be readily ob- 
tained if integral value over arc S tends to the zero. The integral value over 
contour C can be calculated with a help of the Residue Theorem. In some 
cases the evaluating integral will be the real or imaginary part of a contour 
integral. 

Let's consider the function of a complex variable 



g(z)ln(z) 
{zexp(z)) 2 + a 2 
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For this integrand function the conditions of Cauchy's theorem and Jordan's 
lemma are fulfilled. So the integral along the arc of infinite radius is vanishing 
and g(z) have no branch point. Note that all the examples above meet these 
conditions. 

2.1 Integration contour 

Since ln(z) is a multiple- valued function, it is necessary to use the cut plane. 
Usual method (see, for example |§) of dealing with integrals of this type is to 
use a contour large circle C^, center the origin, and radius R; but we must 
cut the plane along the real axis from to oo and also enclose the branch 
point z = in a small circle Co of radius r. The contour is illustrated in 
Fig|IJ. 



Figure 1: Integration contour 
Evidently, we may write down the Cauchy's theorem in a form: 
g(x)\n(x)dx | f g(z)\n(z)dz f g(x)(\n(x) + 2m)dx 



(xexp(x)) 2 + a 2 J (zexp(z)) 2 + a 2 J (xexp(x)) 2 + a 2 

u Coo 00 

r 2 = M £w fa) , (5) 

7 (^exp(2;)) 2 + a 2 ^ 

where is an arc of infinitely great radius, Co - an arc of an infinitesimally 
one. Integrals over these circles are vanishing. Also taking into account the 
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cancellation of integrals evaluated in opposite directions we obtain: 



Thus the evaluating integral is equal to a sum of integrand residues inside 
the contour. The direct calculation of the residue sum in integrand poles leads 
to extremely slow convergence of the partial sums. Nevertheless, we shall see 
below that the convergence series acceleration (Euler's transformation fl2fl ) 
allows to take at most 10 -r- 15 terms into account. 

Following the obvious method one must search the denominator roots and 
substitute the evaluated residues in series (|). 



2.2 Denominator zeros 

The denominator zeros are roots of two equations: 

zexp(z) — ia = 0, and zexp(z) + ia = 0. (7) 

Further we will consider only the first equation of (|7|) as the solution of the 
second one is complex conjugate to the first one. Separating real {x = Rez) 
and imaginary (y = Jmz) parts of equation (|7]) leads to a set (system) of 
two nonlinear equations: 

x cosy — y sin y = 0, ycosy + xsiny — aexp(— x) = 0. (8) 

Simple algebraic transformations allow to reduce this system to: 

xexp(x) . , . 

smy — , x cos y — y sin y = (9) 

a 

As it will be shown below some roots of equation (^) have the positive real 
part (x > 0). Nevertheless, in this case the magnitude xexp(x)/a does not 
exceed unit. For the roots with negative real part (x < 0) at a > 1/e always 
|a;exp(a;)/a| < 1. At a < 1/e there exists a region of negative x values where 
xexp(x)/a < —1. But it's easy to show that no roots of (0) meet this region. 
Let's write down the solution of the first equation (||) in a form: 

y = (-l) Jv arcsin xexp ^ + Nn, N = 0, ±1, ±2, ±3, . . . . (10) 
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Integer number iV separate different branches of function arctg. The y value 
for these branches lie in the limits: 



(2N - 1)~ < y < (2N + 1)|, N = 0, ±1, ±2, ±3, 



After substitution expression ( 10|) into the second equation of the system 
, we obtain: 



/ i \ TV *^ ^-^P /at / ^ 

v — 1) arcsm \- J\ tt = ±a exp{—x / 



\ 



1 - 



xexp(x) 
a 



;i2) 



If A/ = ±2m (m = 0, 1,2, . . .) cosy > and we chose the upper sign 
(plus) in the right-hand side of equation (|12|) . Contrary to this, cosy < if 
A" = ±(2m + 1) and we must chose the lower sign (minus) in the right-hand 
side of equation (|l|). Thus at even A" (N = ±2m) we obtain: 



• xexp{x 
arcsm h i\tt = a exp(— x t 



\ 



xexp(xj 



(13) 



and at odd A" (A" = ±(2m + 1)) respectively: 
xexp(x) 



arcsm ■ 



— Nit = aexp(—x] 



\ 



xexp(x) 



(14) 



It is easy to see, that equation (|T3|) have solutions only at AT > 0, while 
(14) one only at A^ < 0. Therefore, it is convenient to unite (13) and (|i~^) 
and write down: 



. xexp{x 
arcsm V kix — aexp(— x j 



\ 



xexp{x) 



a 



0, k — 0, 1,2, 



(15) 

The point to keep in mind, that cosy > at k = 2m and cosy < 
at k = 2m + 1. At even k = 2m and x < it is not difficult to build an 
approximation for the real and imaginary parts of the root Zk- 



1 + 



In kn I a 
2(/br) 2 _ 



In 



kir 



Vk ~ kir. 



(16) 
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As x < and y > the argument of the complex number is placed in the 
second quarter (less than ir, but greater than tt/2) and is equal to: 

argz k « - + arctg — (17) 
2 m{KTi/a) 

Just note that complex conjugate number argument is placed in the third 
quadrant (x < and y < 0) and: 

37T A/7T 

arg^ fc sa — - arctg (18) 
2 m(k7T/a) 

The comparison of these solutions with exact numerical one will be given 
below in Table 1272 . 




Figure 2: Denominator zeroes in complex plane 

At odd k = 2m + 1 > the argument of the root and complex conjugate 
value are equal respectively: 

37T klT 

arg z k « arctg 



2 & ln(/br/a) 

_ 7T for 
arg2 fc ~ - + arctg . 19 

2 ln(/c7r/a) 
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The comparison of these solutions with exact numerical one will also be given 
in Table £2. 

Hence, evaluation of the complex roots of equation (|7|) is reduced to a real 
part Zk search (transcendental equation (|T5|) solution) and further calculation 
of imaginary part z^ by (|10|) . The graphs of the sum of the first two terms 
of (15) and the third term of the equation (|T5|) are presented on Fig 0. 

Intersection points of these curves marked by circles correspond to the 
roots of equation flT5|). Any of the roots marked on Fig |2|. can be evalu- 
ated numerically with the help of Newton- Rafson method (e.g. [0]). The 
approximations ( |T6| ) are used as an initial guess. Necessary explanation will 
be made in the Program Realization. It is convenient to number these roots 
by corresponding values of k as it is shown on the figure (the numeration is 
used in the program). 



Table 1: Relative error between exact value and approximation 

for different a 



k 


a=l/e 


a=1.0 


a=10 


1 


9.21835- 10" 2 


1.11473- 10" 1 




2 


2.40163- 10~ 2 


2.52883- 10" 2 




3 


1.09204- 10" 2 


1.11533- 10~ 2 




4 


6.21187- 10" 3 


6.27698 • 10~ 3 


6.35174-10" 3 


5 


4.00030 • lO" 3 


4.02305 • 10~ 3 


4.06565 • lO" 3 


6 


2.78844 • 10" 3 


2.79747- lO" 3 


2.82172 • lO" 3 


7 


2.05364-10" 3 


2.05747-10" 3 


2.07187-10" 3 


8 


1.57493- 10" 3 


1.57657-10-3 


1.58551 • lO" 3 


9 


1.24585- 10~ 3 


1.24651 • lO" 3 


1.25229 • lO" 3 


10 


1.01001 ■ 10~ 3 


1.01021 ■ lO" 3 


1.01407- lO" 3 



In case Rez > (it will be if a > \kir\) it is more convenient to find roots 
by bisection because it is difficult to get a necessary initial guess. At the 
same time it is obvious from Fig |2| that at Rez > each subsequent root lies 
between zero and a root found before. 



2.3 Residue Calculation 
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2.3.1 The test integral 

Let's consider the test example of an integral ([3]) when m = 1 and g(r) = 
(1 + r)exp(r). Elementary analytic calculation yields §7r/a. On the other 
hand function F(z) residue in the pole of first order (at the same g(x)) is 
equal: 



resF(z k ) 



[1 + z) exp(z) ln(z) 



d 

dz 



[(z exp(z)) 2 + a 2 



ln(z) 



2zexp(z) 



i ln(z) 



2a 



(20) 



z k 



Zk 



where upper sign (minus) corresponds to the first and lower (plus) to the 
second of the equations (^). 

Let's take the equation (0) in a form: 



[1 + x) exp(x)dx 
(xexp(x)) 2 + a 2 



J2resF(z k ) = — J2M z k) - ln(z k )} 

k=0 Zfl k=0 



i 1 

= 7T ^l iaT &( z k) -iaig(z k )] = — ^][arg(^ fc ) - &rg(z k )]. (21) 
la k=0 la k=0 

For the following evaluation it is important to determine root arguments 
correctly. Here and below in this article the complex value argument will be 
determined by principal value of arctg in order to use the library function 
arctg within the limits of first quadrant. The point to keep in mind is that 
the quadrant of the complex value one have to define "by hand". It allows 
to avoid mistakes connected with possible definitions of library functions of 
negative argument, etc. 

Let's consider first three terms of the series summarized: 

1. If k = the root Zo is placed in the first quadrant and its' argument 
is equal to ipo = arctg \yo/xo\. Evidently, that the complex conjugate 
value has the argument 2tc — ip . Contribution of these two poles in the 
right-hand side of equation (|21|) is: 

arg(^o) - arg(^o) = 2tt - 2ip 



2. The next root z\ (at k — 1) is placed in the third quadrant and its' ar- 
gument is equal to ir + ipi = 7T + arctg \ yi/xi\. Complex conjugate value 
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argument is tt — ipi = tc — arctg \yi/xx\- Corresponding contribution of 
these roots is: 

arg(zi) - arg(^i) = -2^i 

3. Consider one root more at K — 2. Its' argument is equal to 7r — 
{p 2 = n — arctg 1 2/2/^2 1 • For the complex conjugate value we obtain 
7r + <f 2 = tt + arctg \yijx2 1. Thus, we find: 

arg(z 2 ) - arg(^ 2 ) = 2yj 2 
Structure of each term is clear that is why we simply quote the result: 

r (1 + x) exp(x)dx 1 ^ r ... 

/ 7 rw: 2 = o~~ 2J ar S ^ - arg(zfc) = 

J (xexp{x)) z + a A 2a ^ 

1 

[27r - 2(^o - 2(£>i + 2(/? 2 - 2<^ 3 + 2y? 4 — . . .] 



2a 

= - [n - (ip + - if 2 + V?3 - V?4 + • • •)] (22) 
(X 

The sum is found numerically. This series is conventional convergent. 
Partial sums consecutively equal to ~ 0.81 and ~ 2.3 and converge very 
slowly. Nevertheless, the program of convergence acceleration || found the 
result 7r/2 with the relative error that does not exceed 10~ 8 (10 -j- 15 terms 
of the series are taken into account). 

2.3.2 First order pole 

In this section we will consider integral (§) at m = 1 and g(r) = exp(r). The 
residue of integrand function F(z) in the first order pole is evidently equal 
to: 



exp(^) \n(z) 
resF(z k ) = 

— [(zexp(z)) 2 + a 2 ] 



ln(z) 



2z exp(z)(l + z) 



i \n(z) 



2a(l + z) 



*k 

(23) 

As in the previous case the upper sign (minus) corresponds to the first and 
lower (plus) to the second of the equations (0). 
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It is convenient to write down residue expression in a form: 

ln(z k ) ln(z fe ) 



xexp(x)) 2 + a 2 ^ 



2a 



E 

fc=0 



1 + 2 fc 1 + Zfc 



-E 



In \z k \ + iarg(z fc ) In |2 fe | + i arg(2; fe ) 



1 + 



1 + z k 



(24) 



The argument value depends on denominator zeros placement. Obviously, 
there are four possibilities. 

For the complex roots lying in the first quadrant (it will be if k = 
0,2,4..., at a > kit) the first term in ( p4[) is equal to: 

In \z \ + i(p 
l + ^o 

Second term (for the complex conjugate value) yields: 

In \z \ + i(2n - ip ) 
TTio 

Subtracting the terms and separating real and imaginary parts (factor i/2a 
is taken into account) we obtain: 



yo in 




+ 


(1 +x )(tc - ip ) 


a 1 + z 


2 



+ i 



a\l + Zq\ 



(25) 



Pay attention that the angle ipo is determined as the principal value of arctg 
in the limits [0, 7r/2]. 

Similary calculations for the second, third and fourth quadrants give re- 
spectively: 



2/2 


In 


\z 2 \ + (1 + X 2 )(f2 






a 1 + z 2 2 


2/3 


In 


1^3 1 - (1 +^3)^3 






a 1 + 2 3 2 


Zl\ 




(1 + xi)(ir - <fi) 




a 


l + ^i| 2 



+ i 
+ i 
+ i 



ny2 



a\l + z 2 | 2 ' 

7^3 

a|l + z 3 | 2 ' 
a|l + ^i| 2 ' 



(26) 
(27) 
(28) 
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As if structure of each term is already clear then it is not hard to write 
the sum of real residue parts (p5|-|28|) and all the subsequent ones. 



+ 
+ 



exp(x)dx 



(xexp(x)) 2 + a 2 



2/o in 


z 


+ (1 + Xq)(tC - (p ) 




l + z 


2 



J2res(F(z k )) 

k=0 

.:•< , ) 



2/i hi 




| - (1 + Xi)ipx 




1 + Z! 


2 



2/2 In 


Z2 


| + (1 + X 2 )(f 2 


2/4 In 


Z 4 


l + z 2 \ 2 

+ (1 + X 4 )(p4 




1 + Z 4 


2 



2/3 hi 


z 3 


1 - (1 +^3)^3 




l + z 3 


2 



(29) 



Sum of not more than 15 terms of this series give right result for any a > 
0.0006 if acceleration of series convergence is used. The partial sums of the 
imaginary parts also rapidly (~ 10 -v- 15 terms) tends to the zero (~ 10~ 12 ). 



2.4 Approximation at a vanishing 

At a vanishing it is convenient to write down the approximate expression for 
integral (0). Let's set m — 1 in (|3|) and make substitution: 

xexp(x) = au, (30) 

where u is the new independent variable. Left-hand side of equation (|30|) 
may be expanded in powers of the x: 

1 



au ~ x + x 2 + —x 3 + 



Inversion [12] of this series gives: 



x ~ au — (au) 2 + - (au) 3 



(31) 



(32) 



For the present purposes, it is sufficient to retain only one term in right-hand 
side of equation (|3"2|). We will focus on the case discussed in the previous 
section, namely, f(x) = exp(x). Simple integration leads to the following 
result: 



exp(x)dx 



(xexp(x)) 2 + a 2 a(l + a 2 ) 



7T 



+ a ln(a) 



(33) 
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One more term taking into account yields: 

exp(x)dx 
(x exp(x)) 2 + a 2 



7T , . , 3 + 2a 2 , 



v 7 ^- 1 



(34) 



Quality of these approximations is presented in Table |[ 



Table 2: Comparison of approximation and exact value of integral 



a 


Integral value 
(exact) 


relative error for 
approximation (33) 


relative error for 
approximation (|34D 


1Q -07 


1.57079 ■ 10 7 


3.67- 10~ 8 


4.35 • 10" 9 


1Q -06 


1.57078 • 10 6 


3.68- 10~ 7 


4.35 ■ 10~ 8 


1Q -05 


1.57069 • 10 5 


3.67- 10~ 6 


4.35 • 10" 7 


1Q -04 


1.56993 • 10 4 


3.68- 10~ 5 


4.36 • 10~ 6 


6 ■ 10" U4 


2.61115 • 10 3 


2.21 • 10" 4 


2.62 ■ 10~ 5 



As one can see the relative error given by approximation (|34"D is not worth 
than 2.62 • 10~ 5 at a = 6 • 10~ 4 and smaller for other values. If it is necessary, 
the approximation (531) can be improved. 

3 Program Realization 

Program text consists of the principal routine (main) and eight procedures: 

• Root - subroutine for searching real part of denominator zeros. Newton- 
Rafson iteration method used for this purpose. 

• Fun - subroutine for left-hand side equation ([15|) evaluation. 

• Der - subroutine for left-hand side of equation (^) derivative evaluat- 
ing. 

• Approximation - evaluate initial guess to the root according to the (|l6l ) . 
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• Eulsum - program of acceleration of series ( |29| ) convergence. The C 
version of corresponding FORTRAN program || is used. 

• Bisection - recursive version of bisection root finding. 

• Right _bound - auxiliary program. Determine the right bound of root 
region. 

• Sign - auxiliary program for Bisection. Determine sign of variable x. 

The demonstration version of the program is presented below. The first four 

statements declare the necessary functions that are used in the program. 

#include <dos.h> 

#include <stdio.h> 

#include <math.h> 

#include <alloc.h> 

/ / Declaration of functions used below. 

double Fun (double); 

double Der (double); 

double Bisection(double& ,double& ,double& ,double&); 

double Root (double) ; 

double Approximation_x(void) ; 

double Right_bound(void) ; 

void Eulsum(double&, double, int) ; 

int Sign(double) ; 
/ / Defining of the global constants. Parameter A = a can be defined 
// in any other way. EPS - relative accuracy for root search; K - root 

/ / number; N - number of roots taking into account, 
double A = 0.2, EPS = 1.0e-14; 

long unsigned irec = 0, irecmax = 100; 

int K, N = 20; 

/ / Definition of parameters used by procedures Root (x, y, yx, phLk, 
// den, argument) and Bisection (a, c, funa, func, xold, xk). See 
/ / below, 
void main (void) 

{ double x, y, yx, phi_k, den, argument; 
double real_part = 0, imag_part = 0; 
double a, c, funa, func, xold, xk; 

double *real_term = (double*)calloc(N,sizeof (double)) ; 
double *imag_term = (double*)calloc(N,sizeof (double)) ; 
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/ / Two previous statements reserved memory for storaging denomina- 
/ / tor roots. Cycle on root numbers begins below. 

xold = 0.99999*Right_bound() ; 

for(int k = 0; k < N; k++) 

{ K = (k°/ 2) ? -k : k; 

// Variable K (upper case) takes values 0,-1,2 — 3,... while k = 
II 0,1,2,3,.... Statements IF ELSE differs three different cases: 
// Real(zk) = 0, Real(zk) > and Real(zk) < 0. 
if (A == fabs(K*M_PI)) x = 0.0; 
else { if (A > f abs(K*M_PI) ) 

{ a = 0.0; c = xold; 

funa = Fun(a) ; func = Fun(c) ; 
x = Bisection(a, c, funa, func) ; 
xold = 0.99999999 * x; irec =0;} 

else 

{ xk = Approximation_x() ; 

x = Root (xk) ; } 

} 

/ / The following statements calculate argument of z& and store real 
/ / and imaginary parts of the term to be summarized. 

y = pow(-l,k)*asin(x*exp(x)/A) + K*M_PI; 
if(x == 0) argument = M_PI_2; 
else { yx = fabs(y/x); 

phi_k = atan(yx) ; 

argument = (x > 0) ? M_PI - phi_k : phiJs;} 
den = (l+x)*(l+x)+y*y; 

real_term [k] = (0 . 5*y*log(x*x+y*y) + (1+x) * Sign(y) *argument) /den ; 
imag_term [k] =y/ den ; 

> 

/ / End cycle on root numbers. Accelerating the rate of a sequence of 
/ / partial sums performed by the procedure Eulsum. This is the C 
II version of Van Wijngaarden's algorithm (see ||). Then the result 
// is printed. The last two statements free the allocated memory. 

for(int j = 0; j < N; j++) Eulsum(real_part,real_term[j] , j) ; 

for( j = 0; j < N; j++) Eulsum(imag_part , imag_term [j] , j ) ; 
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printf ("Real part =°/,12.51f " ,real_part/A) ; 
printf ("Imaginary part =°/„12.51e \n" , imag_part) ; 
f ree(real_term) ; 
f ree(imag_term) ; 

> 

/ / Initial guess for the Newton- Rafson iteration is chosen with respect 
// to parameter a value. The appointments of other statements are 
evident, 
double Root (double xk) 
{ double xkl ; 

for(int it =0; it < 30; it++) 
{ xkl = xk - Fun(xk)/Der(xk) ; 

if (fabs(xk-xkl) <= fabs(xkl)*EPS) break; 
xk = xkl; } 
return xkl ; } 
/ / No comments, 
double Fun (double x) 

{ double px = x*exp(x)/A, arcsin = asin(px) ; 

return arcsin + abs(K)*M_PI - A*sqrt(l-px*px)/exp(x) ; } 
/ / No comments, 
double Der (double x) 

{ double expa=exp(x), px = A/expa; //px = x*expa/A; 

return (1 + x + px*px) /sqrt (px*px-x*x) ; } 
/ / Calculation of initial guess according to the expressions ([16]) 
double Approximation_x(void) 

{ double mu = abs(K)*M_PI, den = mu*mu, InKPi = log(mu/A) ; 
if(!K) return 0; 

else return -InKPi* (1+0 . 5*lnKPi/den) ; } 

// For details see @. 

void Eulsum(double& sum, double term, int jterm) 
{ double static wksp[28], dum, tmp; 
int static nterm; 

if (jterm == 0) { nterm=0; wksp [0] =term; sum=0 . 5*term; } 
else { tmp = wksp[0]; wksp[0] = term; 
for(int j = 0; j < nterm; j++) 
{ dum = wksp[j+l] ; 



16 



wksp[j+l] = . 5* (wksp [j] +tmp) ; 
tmp = dum; 
} 

wksp [nterm+1] = . 5* (wksp [nterm] +tmp) ; 
if (f abs(wksp [nterm+1] ) <= fabs (wksp [nterm] ) ) 
{ sum += . 5*wksp [nterm+1] ; nterm += 1; } 
else sum += wksp [nterm+1] ; 

} 

> 

/ / Usual bisection algorithm realized by means of the recursion, 
double Bisection(double& left, doublefe right ,double& fun_left, 

doublefe fun_right) 

{ double center = 0.5* (left + right); 
if(++irec < irecmax) 

{if (fabs (left - right) < EPS*fabs (center) ) return center; 
double fun_c = Fun(center); 

if (fabs(fun_left-fun_right)<EPS*fabs(fun_c)) return center; 
if (Sign(fun.left) == Sign(fun.c)) 

{ left = center ; fun_left = fun_c;} 
else { right = center; fun_right = fun_c; } 
center = BisectionQef t , right ,fun_left ,fun_right) ; 

} 

irec — ; 

return center; 

} 

/ / Auxiliary program. Determine the right bound of the root region 
// (see Fig. |). 
double Right_bound(void) 
{ double xn, xnl; 

xn = (A > 1) ? log(A) : A; 
for(int k = 0; k < 20; k++) 
{ xnl = xn + (A*exp(-xn) - xn)/(1.0 + xn) ; 
if (fabs(xnl-xn) <= 1 . 0e-14*f abs(xnl) ) break; 
xn = xnl ; } 
return xnl ; 

} 
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/ / Auxiliary program for Bisection. Determine sign of variable x. 

int Sign (double x) 
{ if(!x) return 0; 

return (x>0) ? 1 : -1; } 

The test result of this programm is presented in Table |3|. 



Table 3: Exact value of the integral and relative error 



a 


Exact 
integral value 


relative error 
for approximation 


1Q -06 


1.57078308845- 10 7 


4.35 • 10~ 8 


1Q -U5 


1.57068696938- 10 6 


4.35 • 10" 7 


1Q -U4 


1.56993298295 ■ 10 5 


4.36 • IO" 6 


1Q -03 


1.56446267294- 10 7 


-2.45 • 10~ n 


1Q -02 


1.53021971263- 10 6 


1.51 ■ 10" 10 


io- 01 


13.7465039696 


1.08 ■ I0~ n 


1 


1.00319691445 


-2.70 • 10~ n 


10 


6.155174431518- 10" 2 


1.05 ■ 10" 10 


100 


3.83403357209 ■ 10" 3 


-1.80 • 10" 10 


1000 


2.62958979905 ■ 10" 4 


-1.04- 10~ 12 



The second column in Table |3| was evaluated by program QUADREC (see 
below) and also checked with the help of MATHEMATICA and MAPLE V. 

4 Conclusion 

Thus it was shown that sufficiently complicated integrals (Q) can be evaluated 
with a given accuracy by means of residue calculations and further evaluation 
of series sum. 

Evidently, that analytical-numeric method like this one, presented in this 
paper, with some restrictions caused by the theory of the function of the 
complex variable may be used for evaluation of a wide class of definite inte- 
grals. 
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A Appendix 

A.l Recursive adaptive quadrature program 

The algorithm consists of two practically independent parts: namely adaptive 
procedure and quadrature formula. 

The adaptive part uses effective recursive algorithm that implements stan- 
dard bisection method. To reach desired relative accuracy of the integration 
the integral estimation over subinterval is compared with the sum of inte- 
grals over two subintervals. If the accuracy is not reached the adaptive part 
is called recursively (calls itself) for both (left and right) subinterval. 

Evaluation of integral sum on each step of bisection is performed by means 
of quadrature formula. The construction of algorithm allows to choose which 
type of quadrature to be used (should be used) throughout the integration. 
Such possibility makes the code to be very flexible and applicable to a wide 
range of integration problems. 

Program realization of the described algorithm is possible when the trans- 
lator that (which) permits recursion is used. Here we propose a C++ version 
of such a recursive adaptive code. The text of the recursive bisection func- 
tion QUADREC (Quadrature used Adaptively and Recursively) is presented 
below: 

void quadrec (TYPE x, TYPE X, TYPE whole) 

{ static int recursion; // Recursive calls counter 

if (++recursion < IP.recmax) // Increase and check recursion level 

{ TYPE section = (x+X)/2; // Dividing the integration interval 

TYPE left=IP. quadrature (x, section) ; // Integration over left subinterval 

TYPE right=IP. quadrature (sect ion, X) ;// Integration over right subinterval 

IP. result += left+right-whole; // Modifying the integral value 

if ( (f abs (lef t+right-whole) >IP . epsilon*f abs (IP . result) ) ) 

/ / Checkup the accuracy 
{ quadrec(x, section, left) ; / / Recursion to the left subinterval 
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quadrec (section, X, right) ;} } 
else IP.rawint++; 
recursion — ; 



// Recursion to the right subinterval 
// Increase raw interval counter 
/ / Decrease recursion level counter 



} 



The form of the chosen comparison rule does not pretend on effectiveness 
rather on simplicity and generality. Really it seems to be be very common 
and does not depend on the integrand as well as quadrature type. At the 
same time the use of this form in some cases result in overestimation of the 
calculated integral consequently leads to more integrand function calls. One 
certainly can get some gains, for instance, definite quadratures with different 
number or/and equidistant points or Guass-Kronrod quadrature etc. 

Global structure IP used in the QUADREC function has to contain the 
following fields: 

struct iip { 

TYPE (*f integrand) (TYPE) ; // Pointers on the integrand and 

TYPE (*quadrature) (TYPE, TYPE) ;// quadrature function 

TYPE epsilon; // Desired relative accuracy 

int recmax; // Maximum number of recursions 

TYPE result; // Result of the integration 

int rawint; // Number of raw subintervals 



The first four fields specify the input information and have to be defined 
before calling the QUADREC function. The next two fields are the returned 
values for the integration result and the number of unprocessed (raw) subin- 
tervals. The static variable RECURSION in the QUADREC function is used 
for controlling current recursion level. If the current level exceeds the spec- 
ified maximum number of recursions the RAWINT field is increased so that 
its returned value will indicate the number of raw subintervals. Here is the 
template for integration of function f(x) over interval [X min , X max ] with the 
use of q(xi,x 2 ) quadrature formula: 

IP . f integrand = f; 
IP . quadrature = q; 
IP. epsilon = le-10; 
IP. recmax = 30; 



} IP; 
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IP. result = IP . quadrature (Xmin,Xmax) ; 
quadrec (Xmin,Xmax, IP. result) ; 



Crude estimation of the integral is evaluated by the quadrature function 
and assigned to IP.result variable. Then it is transferred to QUADREC func- 
tion. Note that the initial estimation of the integral can be set in principal 
to arbitrary value that usually does not alter the result of the integration. 

Further information about QUADREC: test results, using QUADREC in 
MS Fortran source and so on, see [|Il[ . 
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